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Abstract 

We use a hybrid numerical approach to simulate the formation of the Moon from an 
impact-generated disk, consisting of a fluid model for the disk inside the Roche limit 
and an N-body code to describe accretion outside the Roche limit. As the inner disk 
spreads due to a thermally regulated viscosity, material is delivered across the Roche 
limit and accretes into moonlets that are added to the N-body simulation. Contrary to 
an accretion timescale of a few months obtained with prior pure N-body codes, here the 
final stage of the Moon's growth is controlled by the slow spreading of the inner disk, 
resulting in a total lunar accretion timescale of ~ 10^ years. It has been proposed that 
the inner disk may compositionally equilibrate with the Earth through diffusive mixing, 
which offers a potential explanation for the identical oxygen isotope compositions of the 
Earth and Moon. However, the mass fraction of the final Moon that is derived from the 
inner disk is limited by resonant torques between the disk and exterior growing moons. 
For initial disks containing < 2.5 lunar masses (Mcj), we find that a final Moon with 
mass > 0.8Af(j contains < 60% material derived from the inner disk, with this material 
preferentially delivered to the Moon at the end of its accretion. 

Subject headings: Disks, Moon, Planetary formation 
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1. Introduction 

The generally accepted scenario for the formation of the Moon involves the oblique impact of a 
roughly Mars-sized object with the proto-Earth (Hartmann &; Davis 1975; Cameron & Ward 1976; 
Stevenson 1987; Canup 2004). Numerical simulations of such an impact, using primarily Smoothed 
Particle Hydrodynamics (SPH) methods, have shown that the impactor is destroyed (either during 
the impact, or via post-impact tidal disruption), and that a circumterrestrial disk is formed that 
contains up to ~ 2Mj of iron-depleted material (Benz et al. 1986, 1987, 1989; Cameron &: Benz 
1991; Cameron 1997; Canup & Asphaug 2001; Canup 2004, 2008). The sihcate disk is initially a 
mixture of vapor and melt, containing ~ 0(10%) vapor by mass (Canup 2004). Typically about 
~ 20 to 50% of the disk material is predicted to have initial orbits exterior to the Roche limit for 
lunar density material, a^j, with a/j ~ 2.9R^ where is the Earth's radius. 

Prior works have used direct N-body simulations to model the accumulation of the Moon from 
such an impact-generated disk, describing the disk with N = 10^ to 10^ particles that are each 
of order 10^ km in radius (Ida et al. 1997; Kokubo et al. 2000). The N-body models depict an 
extremely rapid disk evolution, with material interior to the Roche limit spreading outward on a 
timescale of order 1 month. Typically a single massive moon accretes in a year or less at an average 
distance of (a) l-3aR (Ida et al. 1997; Kokubo et al. 2000). That a particulate protolunar disk 
would spread so rapidly was anticipated by earlier estimates of Ward &: Cameron (1978), who 
pointed out that a disk of particles containing sufficient mass to produce the Moon would be prone 
to gravitational instability and local clumping. Exterior to the Roche limit, instability-produced 
clumps would form permanent aggregates and seed the growth of the Moon. But interior to the 
Roche limit, such clumps are continually sheared apart by planetary tidal forces, and this process 
generates a large viscosity that drives a ~ lunar mass Roche-interior disk of particles to spread in 
less than a year (Ward & Cameron 1978; Takeda & Ida 2001) (see Section 2.1.1. Both processes 
can be seen directly in the N-body simulations (Kokubo et al. 2000). 

N-body protolunar disk models assume a disk of condensed particles, neglecting the presence 
and creation of vapor as the disk evolves. This is probably a reasonable approximation for material 
orbiting outside the Roche limit. Immediately after the impact, disk material is primarily in 
a condensed state (Canup 2004). Outside the Roche limit, collisions between orbiting particles 
lead to accretional growth. A rough estimate of the heat liberated by accreting the Moon is 
its gravitational binding energy, which implies an energy released per unit mass of the Moon of 
Eb ~ (3/5)GMc/i?c ~ 2 X 10^° erg g~\ where G = 6.67 x 10"® cm^ g~^ s"^ is the gravitational 
constant, and Mc = 7.35 x 10^^ g and = 1738 km are the Moon's mass and radius. Even in 
the limit that all of this energy is retained by the Moon, the expected extent of vapor production 
as the Moon accretes is small because Ei, is much smaller than the latent heat of vaporization of 
silicate, w 2 x 10^^ erg g~^. 

However, there is an inherent inconsistency in describing the Roche-interior disk with an N- 
body particulate model. An approximately lunar mass disk of condensates (solid or melt) interior 



-3- 



to the Roche hmit will be subject to the instability-induced viscosity described above. Such a disk 
spreads so rapidly and the viscously generated heat is so large that the disk would likely substan- 
tially vaporize as it evolves (Thompson & Stevenson 1988), invalidating the model's assumption of 
a particulate disk. A vapor disk would be gravitationally stable, and therefore not be subject to the 
instability-induced viscosity seen in the N-body simulations. Thompson &: Stevenson (1988) were 
the first to recognize this important point, and proposed that the protolunar disk would instead 
evolve in a two-phase, vertically mixed vapor-melt state, with a viscous dissipation rate regulated 
by the rate at which a ~ 2000 K silicate vapor photosphere could cool. This thermally regulated 
viscosity implies a much longer disk spreading timescale of ~ 10^ years (Thompson &: Stevenson 
1988) (see Section 2.1.2). 

Recently Ward (2012) has derived an analytical description for the vertical structure of a 
two-phase protolunar disk inside the Roche limit, including both the vertically well-mixed case 
explored by Thompson &: Stevenson (1988), in which the vapor mass fraction is very low, and a 
new alternative class of solutions in which the vapor initially contains the majority of the disk's 
initial mass. The latter implies a stratified disk structure, in which a portion of the disk's mass 
settles to the mid-plane as melt and undergoes rapid viscous spreading, while the remainder of the 
disk is contained in a gravitationally stable vapor atmosphere. The vapor component of the disk 
requires ~ 10^ years to deplete itself through condensation, and material is ultimately supplied to 
the Roche exterior region over this timescale. Thus both the well-mixed and stratified disk models 
imply a similarly protracted timescale for the inner disk's overall evolution that is ~ 10^ years. 

We here develop a new lunar accretion model that describes the Roche-interior region as a fluid 
disk, while material outside the Roche limit is tracked using direct N-body simulation. The inner 
disk evolves viscously and interacts with outer bodies through resonant torques at the strongest 
Lindblad resonances. Material from the inner disk that viscously spreads beyond the Roche limit 
accretes to form new moonlets that are added to the N-body simulation, while inner disk material 
spreading onto the Earth is removed from the disk. This hybrid construct allows us to model a 
slowly evolving inner disk that spreads in ~ 10^ years as suggested by thermodynamical models 
(Thompson & Stevenson 1988; Ward 2012), while also directly simulating the rapid accretion 
expected among condensed material orbiting outside the Roche limit. 

Our overall objective is a more physically motivated model of the Moon's accretion, and 
improved estimates of its formation timescale and initial orbital position. These quantities are 
related to several outstanding issues, including the Moon's initial thermal state, the potential for 
chemical equilibration between the protolunar disk and the Earth prior to the Moon's accumulation 
(Pahlevan &; Stevenson 2007) , and the initial disk mass and angular momentum required to produce 
a lunar mass Moon. In section 2 we describe in detail our numerical model. In section 3 we use 
our code to reproduce results from pure N-body simulations by Ida et al. (1997). In section 4 we 
perform hybrid model simulations with an inner fluid disk, and study the influence of the disk's 
initial parameters. Results are then discussed in section 5. 
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2. The model 



Our numerical model is built around the symplectic integrator SyMBA (Duncan et al. 1998). 
We have paired it with a simple analytical model for a Roche-interior fluid disk that evolves under 
the influence of viscous spreading and resonant torques due to interactions with orbiting objects 
at their 0*^ order Lindblad resonances. The inner disk mass decreases as material spreads onto 
the planet or outward beyond the Roche limit; mass spreading beyond the Roche limit is assumed 
to accrete into new moonlets that are then added to the N-body code. We assume that the inner 
disk has a uniform surface density a = a{t) and viscosity v = iy{t) with radius. The inner disk's 
evolution is computed by estimating the rate of change of its edges due to the viscous and resonant 
torques. These simplified inner disk treatments are detailed in Appendix A and C. 



2.1. Viscosity model 



We characterize the inner disk by a single, time-dependent viscosity that is a function of the 
disk's surface density a. We envision a silicate disk that is initially two-phase (vapor /melt), and we 
assume that both components co-evolve. We adopt the argument of Thompson &: Stevenson (1988) 
that the inner disk's viscosity will be limited by the rate at which a ~ 2000 K disk photosphere 
can radiatively cool, so long as there is vapor present. Once the disk mass and the associated rate 
of viscous dissipation is low enough that all of the vapor can condense, we assume that the inner 
disk's viscosity will be comparable to that of a purely condensate disk subject to local gravitational 
instabilities (Ward & Cameron 1978). 



2.1.1. Instability induced viscosity 

For a ~ lunar mass, Roche-interior disk composed of melt or solids, local patch instabilities 
strongly increase the collision rate among disk particles through the formation of clumps that are 
continuously sheared apart by planetary tides (Ward & Cameron 1978). By introducing coherent 
particle motions, instabilities modify the transport of angular momentum in the disk and produce 
a characteristic viscosity (Ward & Cameron 1978) 

t^wc ~ — — , (ij 

where a is the disk surface density, and ft is orbital frequency. This process and the resulting rate 
of angular momentum transport are observed in N-body numerical simulations of the protolunar 
disk (Takeda &: Ida 2001) and dense planetary rings (Salo 1995; Daisaka & Ida 1999; Daisaka et al. 
2001). The associated disk spreading timescale for a disk of radial scale r, r^/z^, is then 

rwc = or^o 2 ~ TT y^^"^^- (2) 
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The dissipation rate per unit area is E,y = {9/A)auwc^'^ denotes the time derivative of 
X), implying a total energy per area dissipated as the disk spreads for a time twc of E^twc ~ 
(9/4)(Trf]^. For a uniform surface density disk, this represents a total liberated energy per unit 
disk mass of ~ (9/4)(r$7)^ ~ 5 x 10^^{an/rd) erg g~^. 



2.1.2. Thermally regulated viscosity 

The estimated energy released as an ~ lunar mass disk spreads is thus comparable to the 
latent heat of vaporization for silicate {1^ ~ 2 X 10" erg g~^). If the inner disk spreads in < 1 year 
as implied by equation (2), it is probable that it will retain any viscously dissipated heat, because 
the timescale for the disk to radiatively cool from its surfaces is much longer, of order decades 
(Thompson & Stevenson 1988; Pritchard &; Stevenson 2000). This implies that a substantial fraction 
of a condensate disk would vaporize as it spreads due to an instability-induced viscosity. 

While a condensate disk would be subject to instabilities, a typical silicate vapor disk is 
gravitationally stable. Indeed, its Toomre parameter is (Toomre 1964) 

^ = ^^^%i^J V2000 kJ (lO^gcm-O ' 



where T is the disk's temperature and = ^J^RT/ 7 = 1.4 is the adiabatic index, R is the 
universal gas constant and /u = 30 g mol~^ is the molecular weight. Thus, a protolunar disk 
composed of silicate vapor with a photosphere temperature near the condensation point (T ~ 
2000 K) would have Q > 1, and would not be subject to the instability- induced viscosity. 

In the absence of strong dissipation, a vapor disk could cool and condense, with gravitational 
instabilities then re-developing in the condensed phase. This in turn would heat the disk through 
increased dissipation. This feedback suggests that the dissipation rate in the inner disk will be 
limited by the rate at which the disk can radiatively cool from its surfaces (Thompson & Stevenson 
1988), with 

-^av^^ = 2asBT^, (4) 
and an associated viscosity (Thompson & Stevenson 1988) 

where asB is the Stefan-Boltzmann constant and Tp is the disk's photospheric temperature. The 
corresponding spreading timescale is 
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2.1.3. Model used 

Following the Thompson & Stevenson (1988) disk model in which the liquid and vapor phases 
remain vertically well-mixed, our model assumes that both the vapor and condensed phases vis- 
cously evolve as a single unit. At each time step in our simulation, we compute both the instability- 
induced viscosity i^vKC and the radiation-limited viscosity i'ts (with Tp = 2000 K) at r = r^. If 
i^wc > i^TSi we assume that the disk self-regulates to a radiation-limited viscosity and set v = vts- 
If vwc < ^TS-, the disk can lose energy via radiative cooling at a faster rate than it is generated by 
instabilities. At this point we assume the disk would condense, and therefore set v = I'wc- 

The ratio between the two viscosities is 

VTS 



^2~J3XiU \^2000 Ky' \ajij Vl07gcm-2y " 

For reference, a uniform one lunar mass disk extending from i?^ to aji has a surface density of 
8 X 10^ g cm~^. Our initial disks have vts/^wc < 1 and evolve with a radiation- limited viscosity. 
As the disk spreads and loses mass, its surface density decreases, and since (i'ts/'^wc) cr~^, at 
some point vwc ~ i^TS- This transition occurs for 



4/3 / ^ \ -1/2 



(Ttrans ~ ~ 1-7 X 10^ — ^ — g Cm^^. (8) 



r 



This is equivalent to an inner disk mass of = 1.6 x 10^^ g ~ 0.2 Md for a uniform surface 
density between r = i?^ and = aji. 

We assume that the inner disk maintains a uniform surface density with radius as it viscously 
expands. Viscous expansion leads to mass transfer from the disk onto the planet, and to the 
outward expansion of the disk's outer edge (see Appendix A for details). 



2.2. Spawning of moonlets 

As disk material viscously spreads outward beyond the Roche limit, accretion becomes increas- 
ingly probable and the continuous nature of the disk is disrupted as discrete large objects form (e.g. 
Kokubo et al. 2000). Our model approximates this transition by removing mass from the inner 
disk and adding new moonlets to the N-body simulation once the disk's outer edge expands beyond 
the Roche limit. Once > ur, we compute the characteristic fragment mass that would form 
from local gravitational instability, and assume that since it is at or beyond the Roche limit, it will 
be stable and not be tidally disrupted. This mass, and its corresponding angular momentum, is 
removed from the inner disk and added to the N-body code as a new discrete particle. 

The mass of the fragment that would form from instabilities is (Goldreich &: Ward 1973) 
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where ^ is on the order of, but less than, unity. We set ^ = 0.3. If the disk's outer edge is at 
the Roche hmit {r^ = 2.9R^), then inner disks containing 1.5 and O.OIM^ wiU form fragments 
of ~ 3 X 10~^M(j and 10~^M(j, respectively. This is comparable to the aggregate mass seen in 
the direct "rubble pile" N-body simulations of Kokubo et al. (2000). To improve computational 
efficiency, we set the minimal mass of spawned fragments to ~ 10~^Af(£. Smaller fragments would 
be formed by an inner disk containing < 0.3M{j. As we will later see, this only affects the very 
last stages of a given simulation, so we expect it to be of little influence on the outcome of a given 
simulation. We test the influence of this parameter in Section 5.3. 

At each time step, we check whether the disk's outer edge lies beyond the Roche limit. When 
this occurs, we compute the mass of a spawned fragment per the equation above as a function of a. 
The mass of the spawned moonlet is removed from the inner disk. To conserve angular momentum, 
we first set the new body's semi-major axis to r^, and then we compute the new disk's outer edge 
so that + Lj — = 0, where and are the inner disk's angular momentum before and 
after adding the new body, and Lf is the added fragment's angular momentum. Additional details 
are in Appendix B. 



2.3. Disk-satellite interactions 

We include resonant interactions between the disk and the growing moonlets, which lead to 
a positive torque on the exterior moonlets (whose orbits expand), and a negative torque on the 
inner disk (whose outer edge contracts). Such interactions are important because, e.g., sufficiently 
massive moonlets can initially confine the disk edge within the Roche limit and delay the spawning 
of additional moonlets (Charnoz et al. 2010). 

As a first approximation, we consider only the strongest 0*'^ order inner Lindblad resonances, 
in which the ratio of the mean motion at a location in the disk to that of an exterior moonlet is 
a ratio of integers with (m : m — 1). To compute the resonant torque, we use the formalism of 
Goldreich & Tremaine (1980). 

The total torque Tg exerted by the inner disk on an exterior satellite per unit satellite mass is 
found by summing the torques due to all the 0*'^ order resonances that fall in the disk (see Appendix 
C), 

^""^sGaa,^ C(m), (10) 
where Mg is the satellite's mass, fig = M^/M®, C(m) = ^ 2.55m^(l - 1/m) and is the 

m=2 

highest m for which resonance (m : m — 1) falls in the disk. When m^, < 2, the satellite is 
far enough from the disk that its 2:1 resonance (which is the most distant 0*^' order resonance) 
is no longer in the disk, and in this case Tg = 0. For a satellite orbiting close to the disk, we 
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impose an upper limit on by considering only those resonances that are radially separated from 
the satellite's orbit by a distance greater than the satellite's Hill radius {Rh = as{Ms/3M^)^/^), 
or those resonances for which (1 — < 1 — (Ms/SM^)^^^ . This excludes from the torque 

calculation the approximate region immediately surrounding the satellite's semi-major axis within 
which particles undergo horseshoe orbits. 

To compute the resulting orbital evolution of the satellite we adopt the approach of Papaloizou 
& Larwood (2000). We define an orbital migration timescale, tm, due to the torque on the satellite 
associated with all of its O*'^ order resonances that fall in the disk, tm = Ls/Ts, where is the 
satellite's orbital angular momentum. We then apply an additional acceleration to the satellite, 
given by amig = (v/tm); where v is the satellite's velocity. We include this as an additional "kick" 
of duration r/2 (where r is the timestep) at the beginning and end of each step in the N-body 
code. 

N 

The total torque on the disk due to N orbiting moonlets is = — Tg . We assume that 

s=l 

moonlet torques cause a change in the disk's outer edge r^, with rd\moon < (see Appendix C) 
because external moons remove angular momentum from the disk. 

2.4. Model for Roche exterior particulate disk 

Beyond the Roche limit, we model the protolunar disk material by a collection of individual 
particles, with an initial power-law size distribution N[m)dm oc m~^dm, where N{m) is the number 
of particles with a mass between m and m + dm. In this section we describe our treatments of 
collisions between particles, and the tidal disruption of moonlets scattered close to the Earth. 

2.4- 1- Tidal accretion criteria 

In the default version of SyMBA, all collisions result in inelastic mergers. However this is too 
simplified for objects orbiting near the Roche limit. We modified the code to include tidal accretion 
criteria (Ohtsuki 1993; Canup & Esposito 1995), which near the Roche limit are a function of the 
impact energy, the mass ratio of the colliding objects, and the collision location relative to the Roche 
limit. We use either an "angle-averaged" criterion, that assumes randomly oriented collisions, and 
a "total accretion" criterion, in which collisions are assumed to occur in the radial direction along 
the widest axis of the Hill sphere. These are the same accretion criteria as those used in Ida et al. 
(1997), and in some of the simulations in Kokubo et al. (2000). Additional details can be found in 
Appendix D. 

While an improvement over the assumption of perfect mergers during every collision, our tidal 
accretion criteria are still idealized. They ignore the potential for fragmentation or substantial 
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deformation when calculating whether a given collision results in accretion, and assume that an 
accreted pair merges into a new spherical body. Kokubo et al. (2000) considered three different 
accretion models: the two described above, and a "rubble pile model" , in which individual N-body 
particles are never merged but allowed to form gravitationally bound aggregates of irregular shapes 
that can, e.g., be tidally disrupted when they pass within the Roche limit. Kokubo et al. (2000) 
find similar overall outcomes for all three treatments (e.g., their Figure 1); this is probably because 
the Moon's final position is affected more by its resonant interactions with the inner disk in their 
simulations than the exact position at which it begins to grow, so long as the latter is outside the 
Roche limit. We use the angle- averaged criterion for direct comparison with Ida et al. (1997) in 
our Section 3 simulations, and the total accretion criterion in our hybrid simulations in Section 4. 



2.J^.2. Tidal disruption of moonlets 

Close encounters between particles can lead to some of them being scattered toward the planet 
on high eccentricity orbits, where they may suffer tidal disruption and be effectively absorbed by 
the inner disk. We expect objects accreting in the outer disk to be molten or partially molten (e.g.. 
Section 5.4). An inviscid fluid object on a parabolic orbit will tidally disrupt in a single pass if its 
pericenter distance Q satisfies 



Q < 1.05 ^-^j , (11) 

where Mp is the mass of the disrupting body and po is the density of the orbiting body (Sridhar & 
Tremaine 1992). For the Earth- moon system, this yields Q < 2R^. 

At each time step, we compute the distance of each object to the primary. If this distance 
is smaller than 2i?0, we remove the body from the N-body code and add its mass and angular 
momentum to that of the inner disk. The latter is done by finding the disk's new outer edge 
so that — L'^ + Lc = 0, where L^i and L'^ are the disk's angular momentum before and after 
capture, and is the angular momentum of the captured body. The disk mass after capture is 

= Mfi + rric where rric is the mass of the captured body. This yields 

4 ^5/2 _ ^5/2 ^ ^/5/2 _ ^5/2 

where R = R(q is the disk's inner edge (see Appendix A for details) , and ac and Cc are the captured 
body's semi-major axis and eccentricity. We solve this numerically so that angular momentum is 
conserved to a 10~^ precision. 

To prevent the inner disk's outer edge from expanding too far beyond the Roche limit in a single 
time step due to the tidal disruption of a large object, we implement a tidal stripping mechanism 
for large objects. If a body's mass is greater than 10~^M,5, we remove 20% of its mass at each time 
step once its distance to the primary is r < 2R(q, so that large bodies are entirely disrupted over a 
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few time steps. Since the time step in Symba is 1/20* of the orbital period at 1 Earth radius, large 
bodies are then effectively disrupted over a timescale < 30 minutes, which is much shorter than 
their orbital period. On some runs, this still leads to the outer disk edge temporarily expanding 
to ~ 2.93-R9, but then new bodies are formed from fragmentation (see previous section), and the 
disk outer edge returns to close to the Roche limit in a few tens of orbits. Generally this happens 
late in simulations, when the disk mass is < 10-^ Md, so we believe it does not greatly affect the 
outcome of our simulations. 

We note that the orbits of bodies passing through the inner disk would also be affected by drag 
interaction with inner disk material. For example, an object that encounters a mass comparable 
to its own during a single passage through the inner disk would be captured by the inner disk. 
We neglect this process here, since it depends sensitively on the disk properties (notably its scale 
height and radial surface density profile), which are treated in a simplified fashion by our model. 

3. Tests with pure N-body simulations 

We begin by performing pure N-body simulations, using the angle-average accretion criterion, 
for direct comparison with previous results of Ida et al. (1997), using initial disk parameters given 
in Table 1 of that paper and summarized in our Table 1. 

For runs 1 to 14, the initial disk mass is = 2.44M(j, and the index of the particle size 
distribution, {N{m)dm oc m~^dm) is p = 1.5. Different values are used for the disk's outer edge 
o-maxi the surface density distribution exponent q {cr{a) oc a~'^), and the number of particles N. 
For runs 15 to 19, the disk's outer edge and exponent of the surface density distribution q are held 
constant, while and p are varied. We run each simulation for 5000Tii: where Tk is the orbital 
period at one Earth radius. This timescale is equivalent to that of Ida et al. (1997), who use a 
simulation time of lOOOT^, with being the orbital period at a/j. 

Particles are distributed randomly throughout the disk, with initial eccentricities and inclina- 
tions (in radians) of order O (lO^^) as in Ida et al. (1997). We use their values for the normal 
and tangential coefficients of restitution, with e„ = 0.01 or 0.5 and et = 1. Damping only the 
normal component of the relative velocity of colliding particles can lead to situations in which two 
particles remain in close physical contact but do not actually merge by our accretion criteria, since 
the tangential component of their relative velocity remains unchanged. When we detect such a 
situation, which can cause the simulation's timestep to become prohibitively small, we force the 
merging of the two particles. In practice, in a simulation with ~ 2000 initial particles, this occurs 
between and 2 times. 

The most common outcome of the Ida et al. (1997) simulations was a single large Moon with 
an average semi-major axis of (a) ~ I.Sor. Because inner disk particles spread rapidly and are 
strongly scattered by the forming Moon, all or nearly all of the disk material interior to the Roche 
limit in their simulations was removed in less than a year (through collisions with the Earth or 
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Table 1. N-body simulations parameters. 
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Note. — Input parameters for our pure N-body simulations, adapted from 
Ida et al. (1997). Md, Ld, and Umax are the disk's initial mass, angular momen- 
tum, and outer edge. Ld/Md is the so-called specific angular momentum (in 
units of y/GMf^aa). p and q are the exponents for the initial particle size dis- 
tribution (^N{m) oc m"'') and surface density distribution ((7(a) cx: a"''), and 
A'^ is the initial number of particles in the disk, and et are the normal and 
tangential coefficients of restitution (see text for discussion). Units of mass, 
distance and angular momentum are the present lunar mass M<i, the Roche 
limit for silicates an ~ 2.9i?®, and the angular momentum of the Earth-Moon 
system [Lem = 3.5 x 10*'^ g cm^ 



Table 2. N-body simulations results. 



a 




M a2 


Ma 


M' 


Mpi 


Moo L L' 


Lpi 


Loo 


Run {cir) 


e 


(Mc) (an) 




(Mc) 


(MO 


(Afc) (-Lea/) (iisw) 


(Lem) 


(Lem) 



1 


1 


.54 


0.02 


0.238 


0.60 


0.011 


0.239 


2, 


.160 


0.000 


0.053 


0. 


,054 


0, 


,236 


0.000 


2 





.91 


0.06 


0.120 


1.99 


0.047 


0.286 


2, 


.152 


0.002 


0.021 


0. 


,062 


0, 


,233 


0.001 


3 


1, 


.43 


0.04 


0.367 


0.57 


0.019 


0.369 


2, 


.036 


0.002 


0.079 


0. 


,079 


0, 


,222 


0.001 


4 


1, 


.21 


0.07 


0.162 


2.06 


0.149 


0.322 


2, 


.045 


0.000 


0.032 


0. 


,074 


0, 


,221 


0.000 


5 


1, 


.44 


0.03 


0.424 


6.54 


0.017 


0.442 


1, 


.956 


0.005 


0.092 


0. 


,097 


0, 


,216 


0.002 


6 





.81 


0.16 


0.139 


1.45 


0.081 


0.416 


2, 


.012 


0.012 


0.022 


0. 


,095 


0, 


,222 


0.003 


7 


1 


.23 


0.06 


0.605 


0.61 


0.072 


0.612 


1, 


.693 


0.052 


0.121 


0. 


,123 


0, 


,186 


0.015 


8 


1 


.32 


0.06 


0.734 


0.63 


0.058 


0.737 


1, 


.613 


0.031 


0.152 


0. 


,153 


0, 


,176 


0.010 


9 


1 


.50 


0.14 


0.545 


0.58 


0.179 


0.611 


1, 


.588 


0.061 


0.119 


0, 


,140 


0, 


,176 


0.020 


10 


1 


.19 


0.07 


0.654 


8.93 


0.015 


0.774 


1, 


.591 


0.065 


0.128 


0. 


,165 


0, 


,179 


0.021 


11 


1 


.38 


0.04 


0.827 


3.84 


0.002 


0.830 


1, 


.558 


0.051 


0.175 


0. 


,176 


0, 


,176 


0.017 


12 





.78 


0.28 


0.517 


2.25 


0.505 


1.045 


1, 


.277 


0.119 


0.079 


0. 


,222 


0, 


,144 


0.035 


13 


1 


.57 


0.07 


0.828 


0.96 


0.381 


0.845 


1, 


.096 


0.118 


0.187 


0. 


,192 


0, 


,125 


0.040 


14 


1 


.52 


0.04 


1.315 


0.66 


0.077 


1.315 


0, 


.982 


0.066 


0.292 


0. 


,292 


0, 


,111 


0.023 


15 


1, 


.33 


0.04 


0.537 


0.66 


0.013 


0.551 


1, 


.850 


0.013 


0.111 


0. 


,116 


0, 


,202 


0.004 


16 


1 


.25 


0.02 


0.525 


0.67 


0.016 


0.663 


2, 


.487 


0.063 


0.106 


0. 


,146 


0, 


,275 


0.020 


17 


1 


.40 


0.05 


0.768 


0.56 


0.093 


0.768 


2, 


.321 


0.056 


0.164 


0. 


,164 


0, 


,256 


0.018 


18 


1 


.65 


0.05 


0.496 


0.84 


0.029 


0.504 


1, 


.885 


0.010 


0.115 


0, 


,117 


0, 


,213 


0.004 


19 


1 


.56 


0.06 


0.411 


0.89 


0.018 


0.413 


1, 


.165 


0.000 


0.092 


0. 


,093 


0, 


,129 


0.000 



Note. — a, e M, and L are the semi-major axis, eccentricity, mass, and angular momentum of the largest 
moon at the end of the simulation (t = 5000Tk). aa and Ma are the semi-major axis and mass of the second 
largest body. M' and L' are the mass and angular momentum of the largest body plus all bodies outside 
of its orbit. Mpi and Lpi are the mass and angular momentum of particles that were scattered onto the 
planet. Moo and Loo are the mass and angular momentum of ejected particles. Units of mass, distance and 
angular momentum are the present lunar mass Mc, the Roche limit for silicates ur « 2.9-R®, and the angular 
momentum of the Earth-Moon system (^Lem = 3.5 x lO**^ g cm'^ 
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Moon, or by escape from the system). A primary finding of Ida et al. (1997) (see also Kokubo 
et al. 2000) was that the fraction of the initial disk's mass incorporated into the final Moon is a 
function of the initial specific angular momentum of the disk (L^/M^) and the fraction of the disk 
that escapes during the Moon's accretion (Moo/M^). They used conservation of mass and angular 
momentum to analytically estimate the mass of the largest Moon as: 

M 1.9Ld Moo 

Wd ^ Mdy/GM^an ~ ^'^^ " ^'^ M^' 

where they assumed that the Moon forms at a = l.Sa/j (see also Section 4.3.2). 

Our pure N-body simulations generally reproduce the findings of Ida et al. (1997). In Figure 
1 we plot the ratio of the mass of the largest orbiting body at i = SOOOT;^ to the initial disk mass 
for our simulations, as a function of the disk's initial specific angular momentum. The results 
we obtain are similar to those shown in Figure 5 of Ida et al. (1997). For the most extended 
disks (Runs 13 and 14) we find somewhat larger objects than in Ida et al. (1997), although our 
results for these cases are similar to those obtained by Kokubo et al. (2000) for comparable initial 
Ld/Md values. Analytical estimates from equation 13 with M^o = (solid line) and Mqo = O.OSM^ 
(dashed line) are also plotted on Figure 1. As in Ida et al. (1997) and Kokubo et al. (2000), 
we find that increases for initially more radially extended disks (i.e. for disks with larger 
Ld/Md). The ratio of the Moon's escape velocity to the local escape velocity from the Earth is 
(2GMc/ii(i)^/V(2G'M©/a)i/2 ^ OA{a/aR)^/'^ , so that lunar-sized objects are increasingly effective 
at gravitationally scattering particles into escaping orbits as their orbital radii increase. 

The average semi-major axis, eccentricity, and mass of the final largest moons in our simulations 
are (a) = 1.32a/?, (e) = 0.07, and (M) = 0.54M^, in good agreement with (a) = 1.27aR, (e) < 0.1, 
and {M) = 0.48Mc from Ida et al. (1997). 



4. Simulations with a Roche-interior fluid disk 

We here describe our hybrid simulations that model the Roche-interior disk as a fluid and 
material exterior to the Roche limit with individual particles. For these simulations we adopt the 
total accretion criterion, for which accretion is possible between like-sized objects for a > a/j (see 
Appendix D). 



4.1. Simulation parameters 

Recent impact simulations suggest that the protolunar disk had a mass of 1.5 — 2.1AI<i and a 
specific angular momentum of 0.8 — 1.1, in units of y/GM^OR (Canup et al. 2012). Since we adopt 
a uniform surface density for the Roche-interior disk, the minimum specific angular momentum we 
can achieve, corresponding to a case with only an inner fluid disk extending from 1 to 2.9ii® is 
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Fig. 1. — Ratio of the mass of the largest body at t = 5000Tk to the initial disk mass, as a function 
of the disk's initial specific angular momentum, for pure N-body simulations comparable to those in 
Ida et al. (1997). Squares correspond to Runs with = 0.01 and triangles to those with e„ = 0.5. 
Small symbols are cases where the mass of the second largest body is at least 30% that of the 
largest one, in which cases we plotted the mass of the combined bodies. The solid and dashed lines 
corresponds to equation 13 with M^o = and M^o = O.OSMrf, respectively. 

Ld/Mii w 0.845-^ GM^apt from equation (A3). Simulation parameters are shown in Table 3. We 
consider cases with initial total disk masses = 2, 2.4, 2.5 and 3M^; inner disk masses, Mj„, that 
contain between 50% and 100% of the total disk mass; outer disk edges (global) amax = 2.9, 4, 6, 
7 and 8a r; and exponents for the surface density distribution in the outer disk of g = 1,3 and 5. 
We fix the particle size distribution exponent at p = 1.5, the normal and tangential coefficients of 
restitution e„ = 0.01 and et = l, and the number of particles N = 1500, as those proved to be of 
little influence in pure N-body simulations. An example of the initial setup is plotted in Figure 2a, 
for Run 34. 
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Table 3. Hybrid simulations parameters. 







Li 






Mout 






Run 


yCM^aR) 


{Lem) 


iM<z) 


(Afj) 




q 




1 


0.843 


0.304 


2.00 


2.00 


0.00 


N/A 


2.9 


2 


0.843 


0.365 


2.50 


2.50 


0.00 


N/A 


2.9 


3 


0.955 


0.345 


2.00 


1 


.00 


1.00 


5 


4 


4 


0.960 


0.347 


2.00 


1 


.00 


1.00 


3 


4 


5 


0.965 


0.348 


2.00 


1 


.00 


1.00 


1 


4 


6 


0.955 


0.414 


2.40 


1 


.20 


1.20 


5 


4 


7 


0.960 


0.416 


2.40 


1 


.20 


1.20 


3 


4 


8 


0.965 


0.418 


2.40 


1 


.20 


1.20 


1 


4 


9 


0.899 


0.325 


2.00 


1 


.50 


0.50 


5 


4 


10 


0.901 


0.326 


2.00 


1 


.50 


0.50 


3 


4 


11 


0.904 


0.326 


2.00 


1 


.50 


0.50 


1 


4 


12 


0.899 


0.390 


2.40 


1 


.80 


0.60 


5 


4 


13 


0.901 


0.391 


2.40 


1 


.80 


0.60 


3 


4 


14 


0.904 


0.392 


2.40 


1 


.80 


0.60 


1 


4 


15 


0.888 


0.401 


2.50 


2 


.00 


0.50 


5 


4 


16 


0.890 


0.402 


2.50 


2 


.00 


0.50 


3 


4 


17 


0.892 


0.403 


2.50 


2 


.00 


0.50 


1 


4 


18 


0.880 


0.477 


3.00 


2 


.50 


0.50 


5 


4 


19 


0.882 


0.478 


3.00 


2 


.50 


0.50 


3 


4 


20 


0.884 


0.479 


3.00 


2 


.50 


0.50 


1 


4 


21 


0.986 


0.356 


2.00 


1 


.00 


1.00 


5 


6 


22 


1.009 


0.365 


2.00 


1 


.00 


1.00 


3 


6 


23 


1.036 


0.374 


2.00 


1 


.00 


1.00 


1 


6 


24 


0.986 


0.427 


2.40 


1 


.20 


1.20 


5 


6 


25 


1.009 


0.437 


2.40 


1 


.20 


1.20 


3 


6 


26 


1.036 


0.449 


2.40 


1 


.20 


1.20 


1 


6 


27 


0.914 


0.330 


2.00 


1 


.50 


0.50 


5 


6 


28 


0.926 


0.335 


2.00 


1 


.50 


0.50 


3 


6 


29 


0.940 


0.339 


2.00 


1 


.50 


0.50 


1 


6 


30 


0.914 


0.396 


2.40 


1 


.80 


0.60 


5 


6 


31 


0.926 


0.401 


2.40 


1 


.80 


0.60 


3 


6 


32 


0.940 


0.407 


2.40 


1 


.80 


0.60 


1 


6 


33 


0.900 


0.406 


2.50 


2 


.00 


0.50 


5 


6 


34 


0.909 


0.411 


2.50 


2 


.00 


0.50 


3 


6 


35 


0.920 


0.416 


2.50 


2 


.00 


0.50 


1 


6 


36 


0.890 


0.482 


3.00 


2 


.50 


0.50 


5 


6 


37 


0.898 


0.487 


3.00 


2 


.50 


0.50 


3 


6 


38 


0.907 


0.492 


3.00 


2 


.50 


0.50 


1 


6 


39 


1.068 


0.386 


2.00 


1.00 


1.00 


1 


7 


40 


1.068 


0.463 


2.00 


1.20 


1.20 


1 


7 
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Fig. 2. — Snapshots of the protolunar disk, projected on the R — z plane, at t = 0, 0.03, 1, 30, 
200 and 1000 years, for Run 34 using the hybrid model with a fluid inner disk. The size of circles 
is proportional to the physical size of the corresponding particle. The horizontal thick line is the 
Roche-interior disk. The vertical dashed line is the Roche limit at 2.9i?0. 
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Table 3 — Continued 





Li/Mi 


id 






Mout 




Run 




(Lem) 




(Mc) 


(Ms) 


q {Rs>) 


41 


0.998 


0.361 


2.00 


1.00 


1.00 


5 8 


42 


1.043 


0.377 


2.00 


1.00 


1.00 


3 8 


43 


1.099 


0.397 


2.00 


1.00 


1.00 


1 8 


44 


0.998 


0.433 


2.40 


1.20 


1.20 


5 8 


45 


1.043 


0.452 


2.40 


1.20 


1.20 


3 8 


46 


1.098 


0.476 


2.40 


1.20 


1.20 


1 8 



Note. — Simulation parameters with a Roche-interior fluid disk 
and Roche-exterior individual particles. Md, Ld, and amax are the 
disk's total initial mass, angular momentum, and outer edge. Mi„ 
and Mout are the masses of the fluid disk, and of the solid bodies, 
respectively. Ld/Md is the disk's total specific angular momentum 
(in units of ^GM^au). q is the exponent for the initial surface den- 
sity distribution (o-(a) oc a of the Roche-exterior disk. Units of 
mass, distance and angular momentum are the present lunar mass 
Md, Earth radius R^, and angular momentum of the Earth-Moon 
system (^Lem = 3.5 x 10*^ g cm^ ^~^)- The normal and tangential 
coefficients of restitution e„ and et are set to 0.01 and 1, respectively. 
The particle-size distribution index p is set to 1.5, and the number of 
orbiting particles A'^ is set to 1500. Runs 1 and 2 start with only a 
Roche-interior fluid disk. 
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4.2. Accretion dynamics 



4.2.1. A three-stage accretion 



Figure 3 shows the evolution of the mass of the largest body in Run 34 (solid line), as well 
as the fraction of its mass that consists of material accreted from the Roche-interior disk (dashed 
line). Figure 4 shows the evolution of the number of orbiting bodies. 



Fig. 3. — Mass of the largest body in Run 34 (solid line), and fraction of its mass composed of 
material derived from the Roche-interior disk (dashed line). First Roche-exterior bodies collide 
and accrete until only a few massive bodies remain (1). These bodies confine the inner disk due 
to resonant interactions, and in turn they recede away (2). Eventually, the inner disk viscously 
spreads back out to the Roche limit, and new moonlets are spawned that collide with the Moon 
and complete its growth (3). 

The accretion of the Moon occurs in three consecutive phases: (1) Roche-exterior bodies rapidly 
collide, accrete and scatter one another until only a few massive bodies remain after ~ 1 year (Figure 
2b and c. Figure 3). (2) The inner disk is confined due to resonant interactions with outer bodies, 
which in turn recede away as the inner disk slowly viscously spreads outward. During that time, 
the growth of the Moon is stalled, but the inner disk loses mass on the planet. (3) After ~ 20 years, 
the inner disk spreads back out to the Roche limit, and new moonlets are spawned (Figure 2d). 
These new objects either collide with the Moon to continue its accretion, or get ejected or scattered 
close to the planet where they are absorbed by the inner disk. When the Moon accretes spawned 
moonlets, its semi-major axis tends to decrease slightly, since the specific angular momentum of 
the spawned moonlets is typically smaller than that of the Moon. However interactions between 
the Moon and moonlets that are scattered into the inner disk cause the Moon's semi-major axis to 
increase, as the Moon generally gains angular momentum from the inner scattered bodies (Figure 
2e and f, see also next section). The latter effect dominates the end of the system's evolution. 

Contrary to accretion timescales of less than a year found with pure N-body simulations, here 
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Fig. 4. — Number of orbiting bodies in Run 34. Most bodies initially present in the outer disk 
merge or get scattered in ~ 1 year. After ~ 20 years, the Roche-interior disk has respread back 
out to the Roche limit and starts producing new moonlets as material spreads beyond aR. After 
~ 200 years, the disk produces fragments that are very small, get captured in the 2:1 mean motion 
resonance with the Moon, and are mostly scattered onto the planet (see details in Section 4.2.2). 



the initial confinement of the inner disk by outer bodies and the slow spreading of the Roche-interior 
disk back out to the Roche limit delay the final accretion of the Moon by several hundreds of years. 
We can estimate the minimum mass of an outer object capable of strongly confining the inner disk 
by setting a moon's resonant torque on the disk equal to the disk's viscous torque, assuming that a 
single object confines the inner disk via its 2:1 inner Lindblad resonance. The disk's viscous torque 
at its outer edge reads 

= Snuarjn. (14) 

Assuming that the inner disk contains ~ lM(j, it will evolve with the radiation-limited viscosity 
1/ = UTS (see Section 2.1.3). The confining satellite's torque reads 



3Mffi 



(15) 



Assuming that the confining body's 2:1 resonance lies at the inner disk's outer edge, we have 



rrf = (1 - 1/m) 



2/3 



0.63as, and setting 



requires 



Mr. 



m 



m 



1/3 



1/2 



(16) 



Strongly confining a IMcj disk with rd = an and v = vts requires an outer moon with a mass 
Mf. > 0.07M(j. Thus initially, relatively small moonlets can confine the disk because the radiation- 
limited viscosity is not very strong. Because this viscosity is inversely proportional to the disk's 
surface density, it becomes increasingly difficult to confine the disk as it becomes less massive over 
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time so long as the disk is radiation-limited. However once the disk mass drops to < 0.2M^, the 
viscosity changes to an instability-induced viscosity with = uwc oc cr^, and the disk then becomes 
progressively easier to confine as it dissipates. 



4-2.2. Moonlet-driven orbital migration 

Figure 5 shows the evolution of the mass and semi-major axis of the largest body, and the 
cumulative mass of particles tidally disrupted and absorbed by the inner disk as they are scattered 
close to the planet, for Run 34. While the average semi-major axis of the Moon is ~ 1.3a^ in pure 
N-body simulations, with a Roche-interior fiuid disk this value is increased to ~ 2.15aR (see also 
Table 4). 

8 
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Fig. 5. — Mass and semi-major axis (SMA) of the largest body, and cumulative mass of objects that 
were tidally disrupted and absorbed into the inner disk after being scattered close to the planet. 

The Moon formed here in phase (1) by accretion of the initial outer bodies lies at ~ 4.8i?0. 
During that phase the inner disk's outer edge has been confined within ~ 2.8i2®, so that the Moon 
does not have any resonant interactions with the disk. The latter then slowly viscously spreads 
outward, until it reaches the Roche limit at t ~ 20 years, at which point new moonlets are spawned. 

Initially, the Moon efficiently accretes these moonlets, causing its mass to increase and its 
semi-major axis to decrease slightly (Figure 5, red and green lines). This is due to a change in the 
Moon's angular momentum. Before accreting an object, the latter reads 

L = M^aGM(s (1 - e2), (17) 

where M, a and e are the moon's mass, semi-major axis and eccentricity. After accreting a fragment 
of mass nif, the Moon's angular momentum reads 

L' = {M + mf) yJa'GM(s (1 - e'^), (18) 
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where a' and e' are the post-accretion semi-major axis and eccentricity of the moon. Finahy, the 
fragment's angular momentum is 



~ 4.6i?0, since when it goes inside that distance its 2:1 resonance fahs into the disk, reactivating 
its disk torques and resulting in an outward migration of the Moon. 

In order for a new moonlet to collide with the Moon, their orbits must cross. This can be 
done by an increase in the moonlet's semi-major axis and/or its eccentricity. However, if the latter 
occurs then the object can have a pericenter close enough to the planet that it would be tidally 
disrupted before encountering the Moon. A rapid increase of the moonlet's semi-major axis before 
its eccentricity gets too high is a more favorable scenario for a moonlet to succesfully collide with 
the Moon. 

The inner disk continuously loses mass on the planet and through the Roche limit, so that 
the torque it applies on newly formed objects decreases over time (see Eq. 10). On the other 
hand, the Moon gets more massive over time, making it an even more efficient scatterer. Those 
two effects result in a progressively slower expansion of the semi-major axis of new moonlets, while 
their eccentricity gets excited even more rapidly by the growing Moon. As a result, it becomes 
increasingly difficult for new objects to collide with the Moon before getting tidally disrupted. 
Figure 6 shows the fraction of newly spawned objects that get accreted onto the Moon, tidally 
disrupted, or ejected from the system. It shows that initially, most of the new moonlets collide 
with the Moon. But as time goes by, a larger fraction of objects get scattered toward the planet 
and are tidally disrupted, until this outcome becomes predominant at ~ 120 years. 

A single interaction at conjunction between the Moon and an inner moonlet that is on an ap- 
proximately circular orbit leads to a positive torque on the Moon's orbit. Once the inner moonlet's 
orbit becomes eccentric, subsequent interactions between it and the Moon can lead to a positive or 
negative torque on the Moon's orbit. But if inner moonlets are removed by tidal disruption soon 
after their initial encounters with the Moon, the net torque on the Moon is on average positive, 
which drives an increase in its semi-major axis (Figure 5, green line). When scattering events 
become predominant, the Moon starts migrating outward, at which point it becomes much more 
difficult for objects spawned at the Roche limit to merge directly with the Moon, and thus its 
growth levels off (Figure 5, plateau on the red curve at ~ 120 years). 
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Conservation of angular momentum gives L + Lf = L' . Since the Moon's eccentricity is generally 
of order 10~^ — 10~^, we can set w e'^ ~ 0, which gives 
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Fig. 6. — Fraction of bodies spawned at the Roche hmit that merge with the Moon (black), get 
tidally disrupted (red), or get ejected from the system (green). Initially most bodies merge with 
the Moon, but as the torque from the inner disk decreases due to its decreasing mass, progressively 
more bodies get scattered toward the planet and are tidally disrupted (see text for details). After 
~ 150 years the total of the 3 curves is not equal to 1 because many bodies are trapped in the 2:1 
Moon's resonance, at which point their "fate" has not yet been decided. 

When the Moon's 2:1 resonance is initially located just outside the disk's outer edge at « 3i?®, 
spawned moonlets are not captured into the resonance. The latter requires that the change in a 
moonlet's semi-major axis due to an external torque in one libration period of the resonance be 
much less than the libration width of the resonance (e.g. Dermott et al. 1988). This adiabatic 
condition is violated when the 2:1 is very near the disk edge, because the rate of increase in a 
moonlet's semi-major axis due to disk torques is typically too rapid as it crosses the resonance for 
inner disk masses > 0.1M(j. 

However as the Moon's orbit expands outward due to scattering as described above, its 2:1 
resonance moves away from the disk edge. The disk torque on a moonlet as it crosses the resonance 
is then weaker due to a greater separation between the moonlet and the disk's edge, and capture 
into the 2:1 resonance can occur as the disk is dissipating. For example, in Run 34, newly spawned 
moonlets begin to be trapped in the Moon's 2:1 at about t = 120 years, when the 2:1 has moved 
outward to about 3.1 R®, and the inner disk mass has decreased to ~ 0.5M(j. An inner moonlet 
trapped in the 2:1 resonance continues to receive a positive torque from the inner disk, but because 
it is in resonance with the Moon, the torque causes both the inner body's and the Moon's semi- 
major axes to expand in lock-step. In this way the Moon's orbit is driven outward due to indirect 
resonant interactions with the disk, with the inner moonlets acting as an angular momentum 
relay between the disk and the Moon. Most moonlets captured into the Moon's 2:1 resonance are 
ultimately absorbed by the inner disk and are not accreted by the Moon because the resonance 
prevents close encounters between the moonlets and the Moon, while at the same time increasing 
moonlet eccentricities to high values that lead to close passes by the Earth and tidal disruption. 
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However in cases where multiple moonlets are captured into the 2:1, mutual moonlet interactions 
on occasion scatter objects out of resonance and allow them to be accreted by the Moon. 

Thus due first to the effects of inward scattered moonlets that are lost to the inner disk, and 
then to moonlet capture into the 2:1 resonance with the Moon, the efficiency of accretion onto the 
Moon decreases substantially in the t ~ 120 to t ~ 200 year period, while the Moon's semi-major 
axis increases substantially. This period is also associated with a rapid increase in the total mass 
of particles that are absorbed by the inner disk (Figure 5, black line). 

By t ~ 200 years, the inner disk mass in Run 34 has decreased to ~ 0.2M^, and the 
disk viscosity transitions to the Ward-Cameron viscosity. The viscosity then decreases rapidly as 
the disk dissipates (since i^wc oa a'^), causing the production rate of spawned moonlets to slow 
dramatically (see Figure 7). 

Figure 7 shows the evolution of the Roche-interior disk mass (solid line), of the mass falling 
onto the planet from the Roche-interior disk (dashed lined), and of the mass of new objects accreted 
from the inner disk (dotted line), for Run 34. New objects are accreted at the disk's outer edge 
only after ~ 20 years (Figure 7, dotted line). At that point, the disk has lost ~ 0.04 M(j, that is 
~ 2% of its initial mass (Figure 7, dashed line). At ~ 200 years, the Roche-interior disk is almost 
depleted in mass, and no significant mass is accreted after that point (Figure 7, solid and dotted 
lines) . 
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Fig. 7. — Evolution of the mass of the Roche-interior disk (solid line), of the mass from the inner 
disk fallen onto the planet (dashed line), and of the mass of new objects accreted from the inner 
disk (dotted line), for Run 34. 
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4-2.3. Other possible outcomes 

The aforementioned dynamical steps are present in 28 of the 46 runs, but other outcomes are 
possible. During the first phase of accretion, interactions between initial outer bodies can lead to the 
Moon forming farther away, around ~ 5i?® or beyond. In such cases, the Moon's 2:1 resonance is 
far enough from the disk that newly spawned moonlets are immediately captured in the resonance. 
Subsequent interactions between moonlets can then result in some of them being ejected from the 
resonance and colliding with the Moon, just like in the general mechanism described above, or it 
can lead to the growth of a secondary object inside the resonance. In such cases the Moon can 
be driven out to ~ 8R^, at which point the secondary body's semi-major axis is ~ 4.6i?® where 
it does not interact with the inner disk anymore. The interior body can sometimes get even more 
massive than the body resulting from accretion of the initial outer particles (e.g. Runs 11 and 16). 

For Runs 1 and 2 that do not include an initial outer disk, the outcome is similar to those 
described above. The first body accreted at the Roche limit confines the inner disk inside the 
Roche limit. As the disk viscously spreads outward, the body recoils from the disk and a second 
body is spawned at the Roche limit after ~ 1 year. The disk is once again confined and the second 
body recoils until it finally merges with the first one. This process continues until the Moon gets 
far enough from the disk to start capturing bodies in resonance after several tens of years. At that 
point, interactions between particles leads to the accretion of a second large object that moves in 
resonance with the first one. In Run 1 the resonant configuration remains stable, while in Run 2 
it eventually goes unstable due to interactions with other objects, resulting in a merger of the two 
largest objects. 

4.3. Simulation results 

4-3.1. Properties of the final Moon 

The mass, angular momentum and mass fraction of inner disk material of the largest body from 
each of the hybrid simulations are shown in Table 4. The mass of the largest body at t = 1000 years, 
versus the disk's initial specific angular momentum, is represented in Figure 8. Results show 
a range of outcomes, with an average Moon mass of (M) ~ 0.81 it 0.21 M^, semi-major axis 
(a) = 2.15 ± 0.27a/;, eccentricity (e) = 0.042 it 0.093, and fraction of inner disk material 35 it 30%. 
Accretion efficiency, defined as M/M^, is somewhat lower than in pure N-body simulations, as the 
fraction of the disk accreted varies from ~ 20% to 50%, with the rest being either ejected from 
the system or lost onto the planet. As in Ida et al. (1997) and Kokubo et al. (2000), we find that 
M/Md increases with the initial specific angular momentum of the disk. 
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Table 4. Hybrid simulations results. 
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Morb 


Moo 


Mcap 
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Lorh 


Loo 


Run 


{o-r) 
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(o-r) 


(M<r) 


/2 


(Ms) 


(A^^) 


(Ms) 


(Lem) 


{Lem) 


(Lem) 


1 


2.43 


0.231 


0.370 


100% 


3.88 


0.198 


100% 


0.430 


0.027 


0.330 


0.060 


0.130 


0.011 


2 


2.30 


0.660 


0.016 


100% 


1.39 


0.004 


100% 


0.663 


0.011 


0.483 


0.180 


0.181 


0.004 


3 


1.83 


0.865 


< 10 


3.6% 


1.13 


0.001 


100% 


0.866 


0.017 


0.409 


0.211 


0.211 


0.006 


4 


2.44 


0.714 


0.004 


10.4% 


1.17 


0.001 


100% 


0.716 


0.043 


0.473 


0.201 


0.201 


0.015 


5 


1.87 


0.921 


0.001 


4.4% 


1.16 


0.001 


100% 


0.922 


0.012 


0.472 


0.227 


0.227 


0.004 


6 


1.97 


1.007 


0.001 


5.8% 


1.22 


0.001 


100% 


1.008 


0.026 


0.625 


0.255 


0.255 


0.009 


7 


1.99 


1.066 


0.000 


11.5% 


1.22 


0.002 


100% 


1.068 


0.025 


0.745 


0.271 


0.272 


0.008 


8 


2.07 


0.682 


0.142 


76.4% 


8.46 


0.144 


0.8% 


0.826 


0.035 


0.753 


0.175 


0.244 


0.012 


9 


2.02 


0.703 


0.001 


43.8% 


1.24 


0.002 


100% 


0.705 


0.034 


0.549 


0.180 


0.180 


0.012 


10 


2.05 


0.702 


0.001 


44.5% 


1.29 


0.002 


100% 


0.704 


0.020 


0.558 


0.181 


0.182 


0.007 


11 


2.72 


0.412 


0.026 


15.2% 


1.71 


0.211 


100% 


0.623 


0.027 


0.321 


0.122 


0.171 


0.011 


12 


2.20 


0.798 


0.001 


30.6% 


1.33 


0.003 


100% 


0.801 


0.015 


0.501 


0.213 


0.214 


0.005 


13 


2.18 


0.667 


0.007 


31.2% 


1.36 


0.004 


100% 


0.671 


0.127 


0.563 


0.177 


0.178 


0.041 


14 


1.94 


0.920 


0.002 


38.1% 


1.20 


0.001 


100% 


0.921 


0.022 


0.634 


0.231 


0.231 


0.008 


15 


1.94 


0.925 


0.001 


51% 


1.21 


0.001 


100% 


0.926 


0.030 


0.675 


0.232 


0.232 


0.011 


16 


1.90 


0.397 


0.314 


100% 


3.09 


0.392 


19.8% 


0.789 


0.021 


0.422 


0.093 


0.217 


0.008 


17 


2.02 


0.928 


0.001 


56.3% 


1.22 


0.001 


100% 


0.930 


0.019 


0.689 


0.238 


0.238 


0.007 


18 


2.09 


0.988 


0.002 


52.1% 


1.29 


0.002 


100% 


0.989 


0.047 


0.725 


0.257 


0.258 


0.016 


19 


2.19 


0.998 


0.004 


53.8% 


1.35 


0.003 


100% 


1.002 


0.028 


0.734 


0.266 


0.267 


0.010 


20 


1.84 


0.532 


0.308 


100% 


2.96 


0.354 


27.6% 


0.886 


0.046 


0.703 


0.124 


0.230 


0.017 


21 


2.23 


0.869 


0.001 


8.4% 


1.23 


0.001 


100% 


0.870 


0.022 


0.457 


0.234 


0.234 


0.008 


22 


2.19 


0.865 


0.001 


8.6% 


1.41 


0.003 


100% 


0.868 


0.040 


0.595 


0.231 


0.232 


0.015 


23 


2.20 


0.933 


0.004 


7.2% 


1.34 


0.002 


100% 


0.935 


0.069 


0.546 


0.250 


0.250 


0.024 


24 


2.01 


0.974 


< 10 


6.9% 


1.25 


0.002 


100% 


0.975 


0.118 


0.633 


0.249 


0.249 


0.040 


25 


2.01 


1.054 


0.001 


6.6% 


1.25 


0.002 


100% 


1.056 


0.113 


0.697 


0.270 


0.270 


0.041 


26 


2.02 


1.078 


0.010 


7.6% 


4.54 


0.015 


0% 


1.094 


0.107 


0.709 


0.276 


0.282 


0.038 


27 


2.08 


0.721 


0.004 


37.6% 


1.31 


0.002 


100% 


0.723 


0.021 


0.575 


0.188 


0.188 


0.008 


28 


2.01 


0.770 


0.001 


44.9% 


1.23 


0.002 


100% 


0.772 


0.036 


0.592 


0.197 


0.197 


0.013 


29 


2.78 


0.478 


0.059 


17.5% 


1.73 


0.216 


100% 


0.694 


0.022 


0.220 


0.143 


0.192 


0.009 


30 


1.98 


0.923 


0.001 


45.9% 


1.22 


0.002 


100% 


0.925 


0.044 


0.664 


0.234 


0.234 


0.015 


31 


1.94 


0.919 


0.001 


40.2% 


1.22 


0.001 


100% 


0.920 


0.055 


0.653 


0.231 


0.231 


0.019 


32 


2.96 


0.559 


0.126 


15.6% 


1.86 


0.284 


100% 


0.843 


0.021 


0.290 


0.172 


0.235 


0.008 


33 


1.92 


0.404 


0.316 


100% 


3.04 


0.384 


15.9% 


0.788 


0.037 


0.420 


0.096 


0.216 


0.015 


34 


1.97 


0.951 


0.000 


56.6% 


1.22 


0.001 


100% 


0.952 


0.034 


0.691 


0.241 


0.241 


0.012 


35 


1.99 


0.932 


0.000 


50% 


1.22 


0.002 


100% 


0.933 


0.050 


0.684 


0.237 


0.237 


0.018 


36 


2.86 


0.517 


0.062 


14.3% 


1.75 


0.397 


100% 


0.913 


0.056 


0.441 


0.157 


0.249 


0.019 


37 


2.13 


1.010 


0.002 


56.1% 


1.33 


0.003 


100% 


1.013 


0.077 


0.735 


0.265 


0.266 


0.027 


38 


2.05 


1.107 


0.001 


58.2% 


1.27 


0.002 


100% 


1.109 


0.038 


0.773 


0.286 


0.286 


0.014 


39 


2.05 


0.795 


0.002 


19.3% 


4.27 


0.012 


0.2% 


0.809 


0.215 


0.657 


0.205 


0.210 


0.077 


40 


2.33 


1.067 


0.002 


7% 


1.10 


< 10"^ 


100% 


1.068 


0.097 


0.536 


0.293 


0.294 


0.038 
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Fig. 8. — Ratio of the mass of the largest body at t = 1000 years to the initial disk mass, as a 
function of the disk's initial specific angular momentum. Black, green, red and purple triangles 
correspond to Runs with a total disk mass of 2, 2.4, 2.5 and 3M(j, respectively. Small symbols are 
cases where the mass of the second largest body is at least 20% that of the largest one. In those 
cases we added the mass of the two bodies, since tidal evolution could cause them to merge later 
on (Canup et al. 1999). The black solid and dashed lines correspond to equation (13) with = 
and Moo = O.OSM^, respectively. Blue lines are the analytical estimates from equation (24). 

Figure 9a shows the mass of the largest body at t = 1000 years versus the fraction of its mass 
that is composed of particles accreted from the inner disk. Figure 9b shows the semi-major axis of 
the largest body against its mass. Figure 9c shows the semi-major axis of the largest body against 
its mass fraction of inner disk material. For the less massive initial disks, final moons with a mass 
> O.SMix generally contain less than 20% of their mass in material originating from the inner disk. 
This is because the Roche-interior disk is more strongly confined when a larger object is formed by 
accretion of the initial outer bodies, and while the disk is confined it loses a significant mass onto 
the planet. Starting with initially more massive disks (red and purple points, corresponding to 
total disk masses of 2.5 and 3Mc) increases this fraction to ~ 60%. However lunar- forming impact 
simulations have not generally produced such massive disks for cases in which the impact angular 
momentum is comparable to Lem- 



4-3.2. Analytical estimate 

While analytical estimates from equation (13) are in good agreement with pure N-body simu- 
lations (Figure 1), this is no longer the case for simulations with a Roche- interior fluid disk (Figure 
8, black lines). To derive formula (13), Ida et al. (1997) assumed that the Moon formed at 1.3a/j, 
while in our hybrid simulations (a) ~ 2.15aR. To revise this formula, we redo the calculation in 
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Fig. 9. — a) Fraction of tlie mass of the largest body composed of Roche-interior disk material, 
against mass of the largest body, b) Semi-major axis of the largest body against its mass, c) Semi- 
major axis of the largest body against its mass fraction of inner disk material. Black, green, red 
and purple triangles correspond to Runs with a total disk mass of 2, 2.4, 2.5 and 3M(i, respectively. 
Small symbols are cases where the mass of the second largest body is at least 20% that of the 
largest one. In those cases we added the mass of the two bodies. The two points at 100% inner 
disk material correspond to Runs 1 and 2 that include only a Roche-interior disk initially. 

Ida et al. (1997) and consider conservation of the disk's angular momentum, which gives 

Ld = Mi^GM^ (1 - ef) oi + Mpi^GM^ (l - e^) as + M^^GM^ (l - ej) ag, (21) 

where Mi, ai and ei are the mass, semi-major axis and eccentricity of the IMoon, Mpi = — 
Ml — Moo is the mass scattered onto the planet (solid bodies and through the inner edge of the 
inner disk), with semi-major axis and eccentricity 02 and 62, and M^o is the mass of ejected bodies, 
with semi- major axis and eccentricity 03 and 63. Here we assume that all material initially in the 
disk is either accreted by the Earth, accreted into a IVIoon at o = ai, or scattered onto escaping 
trajectories. Average values from Table 4 give ei PS 0.04 so that (1 - el) ^ 1. Since most of the 
material accreted by the Earth comes from the inner disk (due to tidal disruption, generally no 
particles collide with the Earth), we can set 62 ~ and 02 = -R©. Finally, assuming that the IMoon 
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scatters escaping material on nearly hyperbolic orbits, we can set (1 — 63)03 ~ oi and (1 + 63) ~ 2. 
Equation (21) then becomes 



Ld _ Ml ^ {Md - Ml - Moo)) R® ^ Moo /2ai ^^2) 



Md^GM(saR Md V aR Md \ aR Md \ aR 



Then, using ^yaR/R^ = 1.7 we get 



Ml ^ 1.7 La 1 _ Moo V2^^ - 1 

Md ^ai/R^ - 1 MdVGM^ yJai/R^ - 1 Md T^VoR - 1 ' 



(23) 



For hybrid simulations, assuming the Moon forms at ai ~ 2.15aij, so that ai/R^ ~ 6.2, we get 

^ = 1-14 - 0-67 - 2.3^. (24) 

Md Md^/GM^aR Md ^ ' 

This equation is plotted on Figure 8 as the blue solid and dashed lines, for Moo = and Moo = 
O.OSMrf respectively, which show good agreement with the results from our simulations. 



5. Discussion 
5.1. Summary 

We have developed a new numerical model to study the formation of the Earth's Moon from an 
impact generated disk: an N-body symplectic integrator coupled to a simple model for a fluid Roche- 
interior disk. Our model includes: (1) viscous spreading of the Roche interior disk, using either 
an instability-driven viscosity, or a radiation-limited viscosity, (2) accretion of moonlets when the 
inner disk's outer edge reaches the Roche limit, (3) tidal accretion criteria to treat collisions between 
orbiting bodies, (4) disk-satellite interactions at 0*^ order Lindblad resonances, (5) spawning of new 
moonlets as the inner disk spreads past the Roche limit, and (6) tidal disruption of objects scattered 
close to the planet. 

Our initial setup consists of a fluid disk extending from the Earth's surface to the Roche limit 
at 2.9ii9, and individual particles beyond. We find that the Moon accretes in 3 consecutive phases, 
accreting first from the bodies initially present outside the Roche limit, which confine the inner 
disk within the Roche limit. The inner disk slowly viscously spreads back out to the Roche limit, 
pushing along outer bodies via resonant interactions. After several tens of years, the disk spreads 
beyond the Roche limit, and starts producing new objects that continue the growth of the Moon, 
until the inner disk is depleted in mass after several hundreds of years. For initial disk masses in 
the range impact simulations typically predict, a moon with a mass of 0.6 to 1.1 M(j is produced, 
with a mass fraction of Roche- interior material of 5 to 65%, accreted only during the last stage 
of the Moon's accretion. Increasing the initial total mass of the disk can produced large moons 



- 29 - 



containing up to 60% inner disk material, although it is not clear yet if such disks could be produced 
by appropriate giant impacts. 

Most of our simulations produce a single Moon outside the Roche limit, similar to prior pure N- 
body models (Ida et al. 1997; Kokubo et al. 2000). However there are several key differences. First, 
consideration of a fluid inner disk leads to a lengthening of the Moon's total accretion timescale to 
~ 10^ years, vs. a timescale of < 1 year predicted by pure N-body models. Material that rapidly 
accretes outside the Roche limit resonantly confines the inner disk, which delays the accretion of 
the inner disk material until the disk can viscously spread back out to the Roche limit. For a 
fluid disk with a thermally regulated viscosity (Thompson &: Stevenson 1988), the latter typically 
requires > 50 years. The inner disk material is then preferentially accreted during the last stages 
of the Moon's growth. In contrast, in pure N-body models the viscosity of the inner disk is large 
and the disk spreads very quickly (in < 1 year), so that inner and outer disk material is accreted 
more-or-less simultaneously by the growing Moon. 

The prolonged period of interaction between the fluid inner disk (and moonlets spawned from 
it) and the outer Moon also leads to a substantially larger semi-major axis for the final Moon, with 
(a) ~ 2.15a^ vs. (a) ~ 1.3ajj in the pure N-body simulations. A larger initial semi-major axis for 
the Moon in turn implies a somewhat lower overall accretion efficiency for a given initial disk mass 
and angular momentum, with our hybrid simulations finding that only 20 to 50% of the initial 
total disk mass is ultimately incorporated into the Moon. This suggests that an initial disk mass 
> 2M(j is required to produce a lunar mass Moon, which is somewhat larger than that produced 
to date by most impact simulations (e.g., Canup et al. 2012) that produce a planet-disk system 
whose angular momentum is comparable to that in the current Earth and Moon. 

5.2. Relation to equilibration 

A key constraint on prior impact simulations is the present angular momentum of the Earth- 
Moon system, which constrains the impactor size, the impact angle, and the relative velocity. For 
an impact angular momentum comparable to Lem, the outcome of the impact is a circumterrestrial 
disk composed primarily of impactor material (e.g. Benz et al. 1989; Canup 2004). If the Moon 
accreted from such a disk, it would then have a composition close to that of the impactor. 

The Earth-Moon system shows striking compositional similarities, in particular regarding oxy- 
gen isotopes (Wiechert et al. 2001). However, the distribution of this element in the early solar 
system was very heterogeneous (Clayton 1993). In addition, the scale of radial mixing found in 
terrestrial accretion simulations (Chambers 2001) implies that the impactor would have had a 
substantially different composition from that of the Earth (Pahlevan &: Stevenson 2007), which 
contradicts the observed similarities. 

It has been suggested that mixing could occur between the disk's atmosphere and that of 
the Earth, leading to the equilibration of disk-planet compositions in 10^ — 10^ years (Pahlevan & 
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Stevenson 2007). This is much longer than the accretion timescales predicted by N-body simula- 
tions. However, in our model, the slow spreading of the disk delays the final accretion of the Moon 
by several hundreds of years, which could be compatible with estimated equilibration timescales. 

The 3-step accretion mechanism revealed in our simulations implies that only material accreted 
during the final stage would have been processed through the Roche-interior disk. Earth-like 
material could then naturally end up in the outer parts of the Moon, although mixing in the lunar 
interior would need to be taken into consideration. 

We can however adopt an idealized model for the Moon where initial outer disk bodies accrete 
into a core with radius and material processed in the inner disk piles up later to increase the 
radius to i?2. Noting / the mass fraction of inner disk material, we can express i?i and R2 as 



3M(1-/) 



1/3 



and 



3M 



nl/3 



(25) 



(26) 



For M = 1 Mc, / = 50% and p = 3500 kg m'^, we get Ri ^ 1358 km and R2 ^ 1711 km. 
In the limit that no mixing occurs between the early and late-accreted material, the Earth-like 
material would represent a R2 — Ri ^ 350 km-deep outer layer on the Moon. Whether this would 
be sufficient to explain the identical composition between Earth and the lunar samples is not clear. 



5.3. Model limitations 

5.3.1. Uniform surface density inner disk 

In our model we assume the disk maintains a uniform surface density profile. Numerical 
simulations of the viscous evolution of self-gravitating dense planetary rings show that the disk 
evolves with a density peak inward and lower densities in the outer regions, regardless of the disk's 
initial profile (Salmon et al. 2010). In the instability-driven regime, this would increase (decrease) 
the viscosity close to the planet (at the Roche limit). The opposite would happen in the radiation- 
limited regime (see equations (1) and (5)). A higher viscosity and/or surface density close to the 
Roche limit would decrease the ability of exterior moonlets to confine the inner disk, since the 
balance between the viscous and resonant torque would be more difficult to achieve. As a result, 
more material may be brought viscously through the Roche limit, thus possibly improving the 
accretion efficency, resulting in a larger moon formed for a given disk mass, compared to the slab 
model. However, since the mass necessary to confine the inner disk with the present model is so 
much smaller than a lunar mass, confinement of the inner disk and the associated phase (2) of the 
accretion process seem inevitable when forming a lunar-mass Moon (see Section 4.2.1). Simulation 
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of the radial, as well as temporal, evolution of the inner disk is certainly possible (e.g. Charnoz 
et al. 2010; Salmon et al. 2010), although more computationally intensive, and such modeling is 
planned in our future work. 

Another simplification adopted for the viscous spreading of the fluid disk is the computation 
of the mass fluxes at the disk's inner and outer edges. Both of these fluxes are estimated using the 
viscosity at the Roche limit. However, the viscosity varies with distance as I'wc oc r^/^ or i/ts oc r^. 
Thus we overestimate the rate of mass loss onto the planet, and we expect that future models that 
include the variation of viscosity with distance may increase the disk lifetime. As a consequence, 
we can expect more material to be delivered through the Roche limit, thus increasing the fraction 
of potentially equilibrated material incorporated in the final Moon. 

5.3.2. Co-evolution of liquid and gas phases 

Our inner disk model assumes that both the vapor and condensed phases viscously evolve as a 
single unit. This is motivated by the Thompson & Stevenson (1988) disk model in which the liquid 
and vapor phases remain vertically well-mixed. Recently Ward (2012) has developed a generalized 
description of vertical disk structures appropriate for a two-phase silicate protolunar disk. He 
identifies alternative disk solutions that involve a stratified disk, in which the condensates settle 
to the disk mid-plane and are surrounded by a gravitationally stable vapor atmosphere. The mid- 
plane layer then has a large, instability induced viscosity (per equation (1)), while the atmospheric 
viscosity could be much smaller. The mid-plane layer surface density regulates itself so that the 
energy dissipated matches that which can be radiated from the surface of a ~ 2000i^ vapor disk 
(Ward 2012). 

To describe such a structure will require separate tracking of the condensate and vapor layers, 
which we plan in future work. How might this affect the truncation of the inner disk by the outer 
moon(s)? Initially resonant torques will cause the outer edge of the liquid layer to contract inward 
relative to the outer edge of the gas disk. Gas that lies beyond the outer edge of the liquid layer 
may then condense (because it was the energy supplied by the underlying condensate layer that 
was keeping it in the vapor phase) so that there will be a condensation front that will lie outside the 
liquid layer's outer edge. Once gas has condensed into liquid, the liquid will be subject to resonant 
torques and truncated in a similar manner to that found here. 

5.3.3. Resonances 

For each satellite, we determine which of its resonances fall in the disk and the associated 
torque. The total torque exerted by all orbiting objects onto the disk is then applied at the disk's 
outer edge. As a result, the confinement of the Roche-interior disk may be too efficient. A more 
realistic model that applies torques at the location of each resonance in the disk (as in Charnoz 
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et al. 2010) could increase the ability of inner disk material to spread outward and be accreted by 
the Moon, bringing more potentiahy equihbrated material to the Moon. 

Our model adopts the standard resonant torque expressions appropriate for a cold disk. For 
a hot disk, the positions of the inner Lindblad resonances are shifted inward (and their torques 
correspondingly reduced). A revised torque expression for the (m : m—1) resonance is (Papaloizou 
et al. 2007): 



where ^ is the satellite's gravitational potential, ^ = mh, and h is the disk's aspect ratio, with 
/i ~ 0.1 for the protolunar disk(see Thompson & Stevenson 1988, their Table 1). This expression 
reduces to that given by equation (CI) for ^ ^ 1. 

The high m resonances that are closest to a satellite are the most affected, which will principally 
impact the initial recoil of objects close to the disk's edge. Once a satellite migrates outward away 
from the disk's edge, only its resonances of lower order are in the disk. Of particular importance 
is the 2:1 resonance — which we argue is responsible for the initiation of phase 2 and the nature 
of the Moon's accretion in phase 3 — and the torque due to this resonance is reduced by about 
15% due to thermal effects for h = 0.1. We performed a test simulation using the modified torque 
expression above and found no substantial modifications to the overall accretion history described 
in Section 4.2. 



5.3.4- Size of fragments 

To improve computation efficiency, we set the smallest fragment that can be spawned from 
the inner disk to 1O~^M0, although fragments some two orders of magnitude smaller than this are 
predicted when the inner disk's mass decreases to ~ 10~^M(i. This should not significantly impact 
the outcome of a given simulation, as such small fragments are produced when the disk is almost 
fully depleted in mass, so that no further growth of the moon is expected. To check the influence 
of this parameter, we compare results of simulations using 4 values: no limit, 10~^M®, 10~^M^, 
and IO^^Mq. Table 5 shows the resulting mass, semi-major axis, eccentricity, and fraction of 
inner disk material for the largest and second largest body, at t = 1000 years, for 2 test runs with 
parameters similar to Run 22 (chosen arbitrarily). Although the Moon's predicted semi-major axis 
and eccentricity do increase somewhat as the fragment mass is reduced (due to an increased number 
of particles captured into resonance as their size decreases), results show that the final properties 
of the Moon are not strongly affected by the minimal mass set for the accretion of new moonlets. 
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Table 4 — Continued 





a 


M 






£12 


Ma 




Morb 


Moo 




L 


Lorh 




Run 


(fflii) 




e 


/ 


(ffli?.) 


(Af<r) 


/2 


(Ms) 


(Ms) 


(MO 


(Lem) 


(Lem) 


(Lem) 


41 


2.63 


0.750 


0.050 


13% 


1.64 


0.058 


100% 


0.808 


0.015 


0.297 


0.219 


0.231 


0.006 


42 


2.02 


0.910 


0.007 


12.1% 


4.14 


0.038 


0.3% 


0.949 


0.080 


0.624 


0.233 


0.247 


0.027 


43 


2.00 


0.656 


0.014 


18.2% 


9.02 


0.012 


0% 


0.669 


0.292 


0.828 


0.167 


0.173 


0.106 


44 


2.29 


1.021 


0.001 


8.8% 


1.45 


0.003 


100% 


1.024 


0.060 


0.652 


0.278 


0.279 


0.021 


45 


1.93 


0.897 


0.070 


10.7% 


4.45 


0.221 


0.2% 


1.119 


0.058 


0.847 


0.224 


0.304 


0.020 


46 


2.51 


1.051 


0.002 


9.3% 


1.19 


0.001 


100% 


1.052 


0.122 


0.368 


0.300 


0.300 


0.048 



Note. — a, e, M, f and L are the semi-major axis, eccentricity, mass, mass fraction of inner disk material, and 
angular momentum of the largest Moon at the end of the simulation [t = 1000 years), aa, M2 and /a are the semi- 
major axis, mass, and mass fraction of inner disk material of the second largest body. Morb and Lorb are the mass 
and angular momentum of all orbiting bodies ai t — 1000 years. Moo and Loo are the mass and angular momentum of 
ejected particles. M^ap is the total mass of bodies that where tidally disrupted and captured in the inner disk. Units 
of mass, distance and angular momentum are the present Lunar mass Afj, Roche limit for silicates an ~ 2.97i®, and 
angular momentum of the Earth-Moon system (Lem = 3.5 x 10^^ g cm^ s~^) 



Table 5. Influence of minimal fragment size. 





a 


M 






a2 


M2 






Fragment mass 


(an) 


(Ms) 


e 


/ 


(an) 


(M<r) 


62 






Test Run 1 


Unlimited 


2.82 


0.485 


0.099 


14.60% 


1.76 


0.131 


0.420 


100% 


10-^Me 


2.82 


0.485 


0.100 


14.59% 


1.76 


0.131 


0.420 


100% 


10-*Me 


2.78 


0.485 


0.071 


14.66% 


1.76 


0.134 


0.406 


100% 


10"'' Me 


2.65 


0.487 


0.017 


14.98% 


1.67 


0.141 


0.183 


100% 


Test Run 2 


Unlimited 


2.58 


0.644 


0.043 


26.79% 


1.61 


0.020 


0.557 


100% 


10-^Me 


2.58 


0.644 


0.043 


26.79% 


1.61 


0.020 


0.557 


100% 


IQ-^Me 


2.55 


0.643 


0.005 


26.78% 


1.59 


0.021 


0.283 


100% 


lO-^'Me 


2.56 


0.640 


0.019 


26.39% 


1.59 


0.028 


0.259 


100% 



Note. — a, e, M, and / are the semi-major axis, eccentricity, mass and mass 
fraction of inner disk material of the largest moon at the end of the simulation (t — 
1000 years). 02, 62, M2 and /2 are the semi-major axis, eccentricity, mass and mass 
fraction of inner disk material of the second largest body. Units of mass and distance 
are the present Lunar mass Mc, and the Roche limit for silicates ~ 2.9R^ 
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5.4. Physical state of accreting material 



In the immediate aftermath of the giant impact, ejected material with equivalent circular orbits 
exterior to the Roche limit is predicted to have temperatures ~ 2000K to 5000K (e.g. Canup 2004, 
Fig. 6), and to be predominantly silicate melt (with ~ 0(10%) vapor by mass). If the initial 
outer disk surface density of melt /solid, as, is high enough for gravitational instability, the disk 
will on a rapid, orbital timescale, fragment into clumps with radii R ~ 0(10) km (e.g. Thompson 
& Stevenson 1988). The subsequent timescale for solids in the outer disk to grow through binary 
collisions is 



where p is the bulk density of the solid particles. The timescale for radiative cooling from the 
surfaces of a vertically well-mixed disk is 



where T is the disk temperature, asB is the Stefan-Boltzman constant, and a specific heat Cp ^ 
10^ erg g~^ is assumed. 

Because the time for the disk to cool to temperatures below the solidus is longer than the 
accretion timescale in the outer disk, material initially orbiting outside the Roche limit will accrete 
in a hot, molten state. Pritchard & Stevenson (2000) estimate that protolunar disk material orbiting 
between 2 and 5 Earth radii will take of order 10 years to lose memory of the high temperatures 
produced by the giant impact, and find that individual large objects {R > 100 km) can retain 
temperatures in excess of 1000 K for 10^ years, even given conditions designed to maximize cooling 
(e.g., neglecting the energy of accretion itself). 

Our simulations reveal two stages of accretion: an early, rapid phase in which material initially 
placed outside the Roche limit by the impact accretes in ~ 0.1 years, and a protracted phase in 
which material is delivered to the outer disk on a much longer timescale of ~ 10^ years, as the 
Roche interior disk viscously spreads. During phase 1, accreting objects will be inevitably hot and 
at or above the solidus, while cooling and some solidification might occur during phase 2. 



This work has been funded by NASA's Lunar Advanced Science and Exploration Research 
(LASER) program and the NASA Lunar Science Institute (NLSI). We thank William Ward for 
valuable comments, and for the diffusion model in Appendix A that he developed for Ward &; 
Canup (2000) and Canup & Ward (2000). 
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APPENDIX 
Details on the disk model 



A. Evolution of inner disk mass 

We consider a simple diffusion model for a uniform surface density inner disk by calculating 
the effective changes in its inner and outer edges, R and r^,, under the constraint that a is uniform 
with distance r. This is the same disk model developed by W. R. Ward for Ward & Canup (2000) 
and Canup & Ward (2000). 

The disk spreading timescale is tmsc = ^R^ /v where Ai? = — i? is the disk's width and v 
is its viscosity. Differentiating with respect to time gives 



2 (rj - R) {ri - k) 



(Al) 

where we assume a constant viscosity for simplicity. This assumption is fairly accurate as long as 
integration timesteps remain small, which will be the case here. This yields 

as the rate of expansion of the disk's outer edge. 
The inner disk angular momentum is 

Ld = ^na^GM^ [rf - R^/' 

5/2 _ ^5/2 



= -MdVGM^ , (A3) 

where = air (r^ — i?^) is the disk's mass. Viscous spreading yields no net torque on the disk, 
implying 



which can be rewritten as 



^ ix _ 1) _ 5 (j,2 _ i) j,3/2 



where x = Vd/R- Finally, combining equations (A2) and (A5) gives 



I / A \ 

'^'^'''^^^ " 2R {x-l){l- f{x)) ^ 

2R{x-l){l-fix)y ^'^^> 
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as the rates of change of the disk's inner and outer edges due to viscous spreading. 

We set R = i?®, so that the rate of mass loss from the disk due to infall onto the planet is 



dt 



2-KRRa 



2TrRRMd 



(A8) 



with R from equation (A7). Once expands to the Roche limit, we consider that material that 
diffuses beyond o/j accretes into moonlets that are added to the N-body code (see Section 2.2), and 
this results in an additional loss of mass from the inner disk at a rate 



dt 



2-KrdfdO- 



2TirdrdMd 



1) 



(A9) 



Here = f^d\msc~^ '''''d\moon-> where the first term is from equation (A6) and the second term modifies 
the expansion rate of the outer edge due to satellite torques, per equation (Cll) below. The total 
rate of change in the inner disk mass is then 



dMd 
dt 



dMd 
dt 



+ 



dMd 
dt 



(AlO) 



a,R 



B. Conservation of angular momentum during moonlet spaw^ning 

The angular momentum of the inner disk before fragmentation is given by equation (A3). After 
fragmentation, it becomes 

L'd = -M'dVGM^ _ , (Bl) 
where M'^ = Md — mj. The angular momentum of the fragment is 



Lf=mf^afGM^[l-ej), (B2) 

where af and Cf are its semi-major axis and eccentricity. We set the latter to the ratio of the 
fragment's escape velocity to the local orbital velocity (Lissauer & Stewart 1993) 



if^GM^/a) 



f 



where Rf is the fragment's radius. The fragment's angular momentum then reads 



Lf^m,0aM^(^R,-2^a,y (B4) 

To compute the new disk's outer edge, we numerically solve for L'^ + L f — Ld = so that angular 
momentum is conserved to a 10~^ precision. We then move the new body around its orbit so 
that its actual distance to the new disk's outer edge slightly exceeds its physical radius. Finally, a 
spawned moonlet 's initial inclination is set to half its eccentricity. 



Disk-satellite interactions 



The torque on an exterior moon due to the (m : m—1) inner Lindblad resonance is (Goldreich 
& Tremaine 1980) 



3Qn 



ps 



3nn 



ps 



dr 



2n 



n-n 



ps 



(CI) 



where Q is the orbital frequency in the disk at distance r, ilpg is the pattern speed of the resonance, 
and $m is the m^'^-order Fourier component of the satellite's potential. For O*'^ order inner Lindblad 
resonances, i}ps = i^s where i^s is the satellite's orbital frequency, and the satellite potential can 
be expressed as (Goldreich &; Tremaine 1978) 



——b\/^{a), 



(C2) 



where Ms and are the satellite's mass and semi-major axis, a = r/as = (1 - l/m)2/3 and 6i™^(a) 
is the Laplace coefficient of order 1 /2 defined by 



a) 



cos m 



vr Jq (1 — 2a cos + a^)* 
With il. = mQg/{m — 1), the torque can then be expressed as 



TT a /m — 1 



m 



a, / \ da 



(C3) 



(C4) 



The torque per unit of satellite mass is 



1 



with /is = Mg/M® and Cm = a^^^ I a 
(Goldreich &; Tremaine 1980) 



(C5) 



i m _ T[_ 

Ms~ 3 

+ 2mb^y^ ) . We then use the following approximation 



db 



a- 



(m) 
1/2 



da 



2m 
vr 



A-,|H)+2A-„(| 



2m 

2.51, 

TT 



(C6) 



where Kq and Ki are modified Bessel functions, so that Cm ~ 2.55m^ (1 — 1/m). 

The total torque exerted by the inner disk on an exterior satellite per unit satellite mass is 
found by summing the torques due to all the 0*'* order resonances that fall in the disk, 
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where 



and 



C{m) = c„ 



1 



m=2 

as] 



3/2^ 



(C8) 



(C9) 



where [X\ is the largest integer not greater than X. 

The total torque on the disk due to orbiting moonlets is 



N 



■ Tg . For an inner disk 

s=l 

with a uniform surface density, changing the disk's angular momentum must involve a change in 
its outer edge and/or mass flow across its inner boundary. Because resonances with the outer 
moonlets generally occur in the outer regions of the disk, we assume that moonlet torques cause 
a change in the disk's outer edge r^, with rii\moon < because external moons remove angular 
momentum from the disk, with 



dL„ 



K 3/2 . 



5/2 



dt 5' 

The rate of change of the disk's outer edge due to is then 



^5/2 

2^fdrd 



(CIO) 



rd\ 




(X 



5/2 



1) Ax 



(Cll) 



D. Tidal accretion criteria 

At each time step, we detect which particles are about to collide by checking if their relative 
distance will get smaller than the sum of their radii within the next time step. For each pair of 
colliding particles, we compute the Jacobi energy after the collision (Canup &: Esposito 1995): 

7—1 "^22 ^2 ^^2 ^ ^ /T~\ \ 

Ej = Vi^p - -Xp +2^p- — + 2' 

where Xp, Up and Zp are the Hill coordinates of the impact point, and Vimp is the relative impact 
velocity in units of the Hill velocity Rh^, where = GM^ /og and 
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ao is the local reference radius, nii and m2 are the masses of the colliding particles with physical 
radii ri and r2, rp = (ri + r2)/RH, and e is an effective coefficient of restitution given by 

e = . / W3i, (D3) 

where Cn and €t are the normal and tangential coefficients of restitution, and Vn and vt are the 
normal and tangential components of Vimp- If the post-impact Jacobi energy is < we assume the 
collision will result in a perfect merger (Ohtsuki 1993; Canup & Esposito 1995). 



D.l. Angle-averaged criterion 

When averaging equation (Dl) over all possible impact orientations (radial, vertical, and az- 
imuthal), the post-impact Jacobi energy becomes (Canup & Esposito 1995) 



Ej = le\f^^ - - - + ^. (D4) 



'p 



In the limit that e = (a completely inelastic collision), requiring Ej < for accretion yields 



rp < 0.7. 



D.2. Total accretion criterion 

An alternative is to assume that the particles are aligned in the radial direction, which is the 
widest dimension of the Hill "sphere" and therefore the most favorable for growth. In this case, 
equation (Dl) becomes 

^^ = r"^"mp-^^-lr'p+l m 

In the limit that e = 0, requiring Ej < for accretion then yields < 1. 



D.3. Minimal distance for accretion 

The quantity rp = (ri -|- r2)/RH can be expressed as 



ao V3pcy (1 + ^)1/3 ■ ao(l + ;2)V3' 

where Rc and pc are the radius and bulk density of the central body, p is the material density of the 
colliding particles, or = 2.456i?c(pc/p)^^^ is the Roche limit for material density />, and < /i < 1 
is the mass ratio of colliding particles. Using the above constraints on r^, we can derive a minimum 
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Fig. 10. — Orbital distance (in units of the Roche radius, a^) beyond which two colhding particles 
with mass ratio /x may accrete. The model assumes an inelastic collision between two spherical 
particles. The solid line corresponds to an average over all possible impact orientations, while the 
dashed line considers a purely radial collision (Canup & Esposito 1995). 

distance beyond which two particles of a given mass ratio can accrete, depending on Vimp and e. 
This is represented in Figure 10. 

If the post-impact Jacobi energy of the colliding particles is positive, we assume they rebound 
inelastically and do not merge. The relative velocity of the particles is then modified as 




(D7) 



where v'^ and v[ are the post-impact normal and tangential velocity components. 
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